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ABSTRACT 
Recent  interest  in  military  VTOL/STOL  aircraft  employing  unprepared 
landing  sites  has  led  to  interest  in  the  problem  of  landing  surface  ero- 
sion.  Surface  erosion  is  caused  by  the  aerodynamic  forces  on  ground 
particles  existing  within  the  flow  field  of  an  impinging  jet.   The 
inviscid  flow  field  is  discussed  and  the  viscous  ground  boundary  layer 
is  analyzed  utilizing  both  theory  and  available  experimental  data.   A 
mathematical  model  of  the  process  of  entrainment  of  ground  particles  is 
constructed.   Erosion  rates  in  the  form  of  erosion  profiles  are  pre- 
dicted for  selected  jet  configuration  and  types  of  terrain.   A  criterion 
for  entrainment,  due  to  both  lift  and  drag,  was  found  and  presented  for 
selected  distances  from  the  jet  centerline. 
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NOTATION 

E  erosion  rate  -  ft/sec 

Cr)  coefficient  of  drag 

Cp  coefficient  of  skin  friction 

DN  nozzle  diameter  -  ft. 

d  particle  diameter  -  in. 

F  force  -  lb. 

Lw  lift  force  due  to  the  effect  of  the  ground  -  lb. 

Ls  lift  force  due  to  shear  -  lb. 


m 


mass  -  lb  sec2/ft. 


N  rotational  flow  parameter  =  ^L  /  ,  _ 

NR  Reynolds  number  =  U l  d ^  / 

n  number  of  strips  used  in  strip  analysis 

Ps  static  pressure  on  the  ground  lb/ft2. 

q  shear  velocity  -  ft/sec. 

R  radial  distance  from  jet  centerline  -  ft. 

r  radius  of  particle  -  in. 

T  kinetic  energy  -  ft  lb. 

U  velocity  parallel  to  the  ground  -  ft/sec. 

V  velocity  perpendicular  to  the  ground  -  ft/ sec. 

x  coordinate  parallel  to  the  ground 

y  coordinate  perpendicular  to  the  ground 

y,  vertical  distance  at  which  U  =  Urnag.        in  turbulent  boundary  layer 
^  z. 

Z  height  of  nozzle  above  the  ground  -  ft. 

\  density  -  lb  sec^/ft  . 

/M  viscosity  -  lb  sec/ft2. 


©  angle  measured  from  center  of  particle 

6J>o  vorticity 

O  boundary  layer  thickness 

Subscripts 

crit.  value  at  boundary  layer  transition 

J  indicates  value  associated  with  jet 

L  value  at  local  point  on  the  particle 

a  property  of  air 

I  value  in  inviscid  flow 

p  value  associated  with  a  particle 

c  measurements  made  on  cylinder  in  strip  analysis 

s  static 

A  ambient 

n  measurement  made  on  a  strip 

w  value  associated  with  effect  of  wall  constraint 

S  value  associated  with  shear  effect 


vi 


1.    Introduction. 

The  increasing  interest  in,  and  use  of,  VTOL  aircraft  and  other 
vehicles  with  vertically  directed  fans  and  jets,  particularly  in  mili- 
tary operations,  has  brought  about  concern  for  the  deterioration  of 
unprepared  or  soft  landing  sites.   The  downwash  impingement  caused  by 
these  aircraft  results  in  a  radially  developed  viscous  flow  field  on  the 
ground,  producing  forces  adequate  to  cause  entrainment  of  the  ground 
particles,   The  reasons  for  concern  with  this  problem  are  numerous.   In 
addition  to  possible  damage  to  the  aircraft  itself  and  the  erosion  of 
the  landing  surface  there  are  the  accompanying  hazards  to  ground  per- 
sonnel, damage  to  ground  support  equipment,  poor  visibility  afforded  the 
pilot,  and  the  violation  of  the  military  concept  of  concealment.   The 
present  study  was  designed  to  construct  a  mathematical  model  of  the 
ground  erosion  process  to  allow  prediction  of  the  initial  erosion  rates 
for  various  types  of  terrain  for  various  nozzle  velocities  and  geome- 
tries.  Previous  work  \\,    14,  15,  22,  23,  25J  on  the  problem  of  erosion 
due  to  downwash  impingement  has  been  primarily  an  experimental  examina- 
tion of  the  problem  with  small  scale  jets. 

The  critical  mechanism  in  the  erosion  process  is  the  interaction  of 
the  ground  boundary  layer  and  the  terrain  particles.   An  analysis  of 
this  mechanism  requires  an  understanding  of  the  viscous  mixing  problem 
of  a  finite  free  jet,  with  the  inviscid  flow  field  extending  radially 
from  the  jet  center  line  as  a  result  of  the  presence  of  the  ground  plane, 
An  analysis  of  the  phenomenon  which  produces  entrainment  of  the  ground 
particles  is  basically  a  study  of  the  aerodynamic  forces  which  exist  in 
the  ground  boundary  layer,  as  they  are  affected  both  by  the  decaying  in- 
viscid  flow  field  above  the  boundary  layer  and  by  the  constraint  of  the 


ground  surface.  In  analyzing  the  boundary  layer  forces  available  theory 
was  utilized  whenever  possible,  but  since  the  radially  developed  bound- 
ary layer  runs  the  gamut  of  stagnation  point  flow  to  turbulent  flow  and 
then  complete  decay,  there  are  areas  for  which  no  theory  or  exact  solu- 
tion exists.  To  fill  these  holes,  experimental  data  were  used.  Where 
exact  solutions  were  used  for  the  velocity  profiles  within  the  boundary 
layer,  verification  was  made  with  experimental  results  if  data  existed. 

The  analysis  was  restricted  to  single  jets  of  uniform  dynamic  pres- 
sure loading  and  circular  cross-section  under  "no  wind"  conditions  shown 
diagrammatical ly  in  Fig.  1.   An  axi-symmetric  boundary  layer  was  assumed, 
acting  over  a  terrain  of  uniform  particle  size  and  density.   Table  I  in- 
dicates the  various  combinations  of  nozzle  characteristics  and  radial 
stations  that  were  analyzed.   The  prediction  of  erosion  rates  is  neces- 
sarily restricted  to  incipient  motion,  since  the  secondary  effect  of 
random  collision  is  not  readily  treated  by  analysis  and  because  the  con- 
stantly changing  shape  of  the  ground  surface  is  not  accounted  for  as  a 
continuing  process.   Although  this  work  is  restricted  to  consideration 
of  uniform  jets  only,  previous  experiments  by  Morse  \_23j  indicate  that 
the  dynamic  pressure  distribution  due  to  downwash  from  rotor  blades  and 
ducted  fans  is  similar  to  that  for  a  uniform  jet.   The  uniformity  of 
the  jet  or  downwash  appears  to  be  a  critical  factor  only  for  ratios  of 
jet  altitude  to  diameter,  Z/D^j,  less  than  unity. 


2.    History  of  the  Problem. 

The  classic  problem  of  an  impinging  jet  on  a  normal  surface  is  not 
new,  but  its  application  to  the  field  of  VTOL  aircraft,  of  either  rotor 
or  jet  type,  has  been  generated  only  within  the  past  several  years.   The 
problem  of  surface  erosion  resulting  from  impinging  downwash  first  re- 
ceived attention  from  Kuhn  [l4].   Rutin' s  work  was  experimental  employing 
small  scale  test  equipment  to  make  dynamic  pressure  distribution  surveys. 
An  experiment  was  conducted  by  fixing  the  nozzle  dynamic  pressure  and 
raising  or  lowering  it  over  terrain  of  known  condition.   Testing  was 
done  with  various  types  of  sod,  dry  and  wet  dirt,  dust,  and  sand.   Simi- 
lar experimental  tests  have  been  conducted  more  recently  by  Morse  \_23J . 

The  most  recent  work  done  on  this  problem  is  that  of  Brady  and  Lud- 
wig  [26] .   This  work  consists  of  extensive  experiments  to  investigate 
the  characteristics  of  the  impinging  uniform  jet,  including  both  the 
inviscid  flow  and  the  ground  boundary  layer.  Measurements  of  boundary 
layer  velocity  profiles  were  compared  with  theory  and  showed  agreement 
to  a  limited  degree,,   Some  of  these  experimental  results  have  been  used 
in  support  of  the  present  work.   Vidal  [28]  has  given  an  insight  into 
the  aerodynamic  forces  existing  within  the  boundary  layer. 


3.    The  Flow  Field. 

In  order  to  consider  the  mechanism  of  particle  entrainment  and  ground 
erosion,  a  thorough  knowledge  of  the  flow  field  of  the  impinging  jet  is 
required.   This  flow  field  can  effectively  be  separated  into  three  re- 
gions.  The  first  of  these  regions  exists  regardless  of  the  presence  of 
the  ground  surface,  and  is  the  region  of  viscous  decay  immediately  exter- 
nal  to  the  jet  nozzle.   Viscous  decay  of  a  free  jet  is  a  classical  prob- 
lem under  uniform  conditions.   This  region  is  characterized  by  the  flar- 
ing out  of  the  jet  as  the  jet  boundary  mixes  viscously  with  the  stationary 
field  external  to  it.   Introducing  a  ground  surface  into  the  flow  field, 
normal  to  the  jet,  does  not  invalidate  the  quantitative  solution  of  free 
jet  decay,  but  restricts  its  extent  of  usefulness.   Experimentation  has 
shown  that  the   theory  of  viscous  decay  of  a  free  jet  is  adequate  for 
the  impinging  jet  problem  for  nozzle  heights  of  greater  than  about  eight 
jet  diameters.   Since  in  this  work,  the  interest  is  in  nozzle  heights 
of  Z./o   <_  <Q     ,    an  additional  region  of  viscous  mixing  with  the  station- 
ary boundary  must  be  studied  and  determined,  beginning  with  the  point  at 
which  the  presence  of  the  ground  surface  has  altered  the  free  jet  charac- 
teristics.  Further  study  must  then  be  made  of  the  viscous  decay  of  the 
jet  after  impingement  on  the  ground  plane.   To  date  no  satisfactory 
theoretical  method  exists  for  dealing  with  the  entire  flow  region  in 
three  dimensions.   This  problem  may  be  partially  overcome,  however,  by 
including  the  viscous  decay  region  in  describing  the  inviscid  flow  field. 

Knowledge  of  the  inviscid  flow  field  is  necessitated  by  the  require- 
ment that  the  inviscid  velocity  existing  at  the  upper  edge  of  the  ground 
boundary  layer  be  known.   The  difficulty  in  finding  a  theoretical  solu- 
tion for  this  region  lies  in  the  fact  that  although  the  boundary  conditions 


are  well  defined,  the  location  of  the  free  streamline  is  not  known.  Three 
dimensional  theoretical  solutions  of  this  problem  are  limited  to  approxi- 
mate methods.   One  of  these,  by  LeClerc  ^JLOQ  >  uses  an  electrolytic  analog 
to  determine  the  shape  of  the  free  streamline  boundary.   In  this  method 
the  solution  of  the  inviscid  flow  region  was  accomplished  by  relaxation 
techniques,  predicting  velocity  and  pressure  distribution.   Such  a 
method  of  solution  could  be  made  to  approach  the  exact  solution  to  any 
desired  degree  of  accuracy. 

In  the  present  work  an  exact  solution  of  the  inviscid  flow  field 
was  not  required.   It  was  assumed  adequate  to  use  existing  experimental 
data  which  compared  favorably  with  the  limited  theoretical  analyses  avail- 
able.  In  most  instances  it  was  found  that  the  theoretical  solutions  were 
not  in  good  agreement  with  the  experimental  data  for  (VoN<zJ-.   The  approx- 
imate methods  of  predicting  the  characteristics  of  the  inviscid  flow  re- 
gion by  use  of  an  equivalent  inviscid  jet  of  greater  diameter  and  decreas- 
ing dynamic  pressure,  or  by  replacing  the  jet  with  a  cylindrical  vortex 
sheet  of  constant  radius  extending  to  the  ground,  have  proven  unsatis- 
factory when  compared  with  experimental  data.   In  view  of  this,  available 
experimental  data  were  used  in  this  work  as  a  solution  to  the  inviscid 
flow  field.   Fig.  2  shows  the  result  of  empirical  equations  fitted  to 
these  data. 

The  third  flow  region  is  that  of  the  ground  boundary  layer.   A  thor- 
ough knowledge  of  this  region  was  most  important  to  this  work  as  the  ground 
particles  are  predominately  immersed  in  the  boundary  layer,  and  the  mechan- 
ics of  particle  movement  originate  with  the  aerodynamic  forces  resulting 
from  the  characteristics  of  the  boundary  layer. 


4.    The  Ground  Boundary  Layer. 

The  key  to  the  entire  analysis  of  ground  erosion  is  the  mathematical 
model  used  for  the  boundary  layer.   Sound  theoretical  analysis  of  the 
boundary  layer  was  used  in  this  study  when  possible.    Experimental 
results  were  used  whenever  the  theory  was  non-existent  or  was  in  large 
disagreement  with  these  results.   For  effective  analysis  the  boundary 
layer  may  be  divided  into  four  regions.   The  first  is  the  stagnation 
area  boundary  layer  existing  adjacent  to  the  jet  center  line  and  extend- 
ing to  R/D  =z  .5  or  to  a  station  directly  under  the  original  free  jet 
boundary.   This  area  is  laminar  for  the  conditions  investigated.   The 
solution  of  this  region  utilized  the  Himenez  solution  t_18 J  for  the  uniform 
impinging  jet  on  a  perpendicular  wall.   Experimental  results  were  not 
available  for  this  region,  so  that  the  validity  of  this  solution  is  not 
verified  but  was  assumed.   The  Himenez  solution  in  three  dimensional 
flow,  developed  by  Homann  [l8j  is  an  exact  solution  utilizing  the  Navier- 
Stokes  equations  for  axisyrametric  flows.   The  solution  was  developed 
explicityly  for  stagnation  conditions  and  the  velocity  profiles  result- 
ing from  this  have  been  employed  to  R^  •=  0.5.   The  nondimensional 
velocity  profile  is  shown  in  Fig.  3. 

The  second  region  is  the  laminar  boundary  layer  extending  from  the 
edge  of  the  original  jet  boundary  to  the  radial  station  at  which  transi- 
tion begins.   A  complete  method  of  solution  for  this  laminar  region,  by 
Smith  [l9]  ,  is  well  verified  by  experiment  for  ^-/q^  =  .5  [26],  but 
gives  little  insight  as  to  the  limit  of  the  laminar  region  and  the  be- 
ginning of  transition.   This  point  has  been  approximated  by  Brady  and 
Ludwig  [.26]  from  considerations  of  neutral  stability  of  the  laminar  bound- 
ary layer.   A  Reynolds  Number  at  which  the  boundary  layer  becomes  neutrally 


stable  can  be  determined  [18J  and  may  be  converted  to  a  jet  nozzle 
Reynolds  Number  as  a  function  of  V^/d^i")  crit:  '   0n  tn*-s  basis  the 
boundary  layer  becomes  unstable  at  'Yo^^^  1*0  for  the  jet  velocities 
and  nozzle  diameters  analyzed  in  this  work.   In  view  of  this,  and  the 
fact  that  the  invisced  velocity  generally  reaches  a  peak  in  the  vicinity 
of  ^y^  =  1.0  (Fig.  2),  the  beginning  of  transition  was  taken  to  be 
R/Drs|  =   1.0. 

The  transition  region  is  even  more  difficult  to  define.   Since 
fully  turbulent  flow  exists  theoretically  when  the  pressure  gradient  on 
the  ground  surface  is  essentially  zero,  it  was  assumed  that  the  transi- 
tion region  extended  to  the  station  where  the  pressure  gradient  could  be 
taken  to  be  negligible.   Fig.  4  shows  the  static  pressure  distribution 
near  the  ground  taken  from  an  electrolytic  analogy  of  the  entire  flow 
field  ^10]|.   The  station  at  which  the  gradient,  &\  r~^j-^J  /  \  y>      } 

was  taken  to  be  zero  was  *Vdn  =  2.0.   For  the  laminar  and  transition 
regions  of  the  boundary  layer  experimental  results  were  used  to  form  the 
velocity  profiles.   This  data,  collected  by  Brady  and  Ludwig  [26]  for  a 
limited  range  of  jet  velocities  indicate  that  the  non-dimensional 
velocity  profile  was  almost  completely  independent  of  the  jet  velocity. 
The  measurements  were  taken  at  nozzle  heights  in  the  range  ^/d^"  0*5 
to  4.0  at  four  locations  and  for  the  radial  range  of  ^/oN\=«5  to  1.33 
at  six  different  stations.   These  twenty- four  configurations  are  assumed 
to  be  fairly  representative  of  the  laminar  and  transition  regions  and 
are  shown  in  Figs.  5a  through  5h.   For  purposes  of  computation,  curves 
were  fitted  to  each  group  of  data. 

The  fourth  region  in  the  boundary  layer  is  the  turbulent  regime. 
It  was  assumed  that  this  region  begins  at  ^/c.^^  2  as  this  is  the 
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station  at  which  the  pressure  gradient  is  assumed  negligible.   A  limited 
amount  of  work  has  been  done  on  theoretical  solutions  of  turbulent  bound- 
ary layers.   One  of  these  that  appears  applicable  to  the  impinging  jet 
is  that  of  Glauert  [_7j.   This  single  solution  for  velocity  profiles  in 
the  turbulent  boundary  layer,  assumed  valid  to   /oN-  4.  which  is  the 
radial  limit  of  the  analysis,  has  been  verified  by  experimental  results 
at  R/oN  -  1.33  [26].   The  uncertainty  as  to  the  characteristics  of  flow 
decay  limit  further  extension  of  the  analysis. 

Glauert  set  up  the  boundary  layer  equations  and  found  a  similarity 
solution  in  which  the  form  of  the  velocity  distribution  across  the  jet 
does  not  vary  along  its  length.   The  solution  is  dependent  upon  the 
assumption  of  an  effective  eddy  viscosity  as  required  to  satisfy  the  law 
due  to  Blasius  for  flow  in  a  pipe  as  well  as  Prandtl's  hypothesis  for 
free  turbulent  flow.   The  first  solution  is  considered  to  be  valid  near 
the  ground  surface,  and  the  second  to  be  valid  above  some  determined  dis- 
tance from  the  ground  surface.   Fig.  6  shows  these  results  in  the  form 
of  a  velocity  profile,  for  the  radial  distances  and  nozzle  heights  con- 
sidered in  the  present  work.   It  was  found  that  the  shape  of  the  velocity 
profile  was  essentially  independent  of  both  the  radial  distance  and  the 
nozzle  height,  but  is  dependent  upon  the  nozzle  Reynolds  Number. 

The  turbulent  region  is  characterized  by  a  velocity  profile  which 
approaches  a  maximum,  U  r^o^x  ,  with  increasing  y  and  then  decreases  with 
further  increase  in  y.   As  is  shown  in  Fig.  1  it  is  expected  that  the 
turbulent  boundary  layer  will  continue  to  grow  as  long  as  there  exists 
an  inviscid  upper  boundary.   Since  very  little  attention  has  been  given 
to  this  phenomenon  of  flow  decay,  and  because  of  the  nature  of  the  tur- 
bulent solution,  a  fictitious  boundary  layer  thickness  is  assumed  for 
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this  region.   This  is  defined  to  be  the  thickness  that  exists  at 

\J  -  Urr\ooc  *   *fc  was  f°unc*  that  the  nature  of  this  assumption  does  not 

affect  the  overall  results  appreciably. 


5.   Analysis  of  Aerodynamic  Forces  Within  the  Boundary  Layer. 

The  ground  boundary  layer  is  a  non-uniform  flow  that  is  character- 
ized by  a  large  velocity  gradient  or  shear  layer  extending  over  a  large 
portion  of  the  thickness.   In  addition,  the  ground  surface  or  wall  bound- 
ary at  the  base  of  the  flow  forms  a  constraint  on  the  streamlines  as 
they  attempt  to  pass  beneath  an  obstacle  bedded  in  the  ground  surface. 

There  are  three  distinct  forces  that  exist  within  the  boundary 
layer.   These  are:   1)  Drag  due  to  the  horizontal  velocity  component  of 
the  flow;  2)  lift  due  to  the  proximity  of  the  wall;  and  3)  lift  due  to 
the  velocity  gradient  in  the  flow.   There  is  no  theory  available  which 
deals  with  this  problem  in  general.   Therefore  it  was  necessary  to  piece 
together  an  approximate  theory  for  each  of  the  various  forces. 
Drag  Force. 

The  drag  force  was  the  easiest  to  account  for.   To  simplify  the 
analysis  in  a  non-uniform  flow,  a  strip  analysis  was  made  in  which  the 
spherical  particle  was  divided  into  a  number  of  layers  and  it  was  assumed 
that  a  uniform  flow  acted  upon  each  layer.   Knowing  the  velocity  of  the 
flow  on  each  layer,  from  the  velocity  profile,  the  force  on  each  layer 
and  the  summation  of  horizontal  forces  were  developed.   The  calculation 
of  the  total  drag  force  employed  the  drag  coefficient  of  cylindrical 
bodies  measured  by  Hoerner  £sT): 


N\r^    <   i  oo 

-      30,/ 


<-*     -    3%>^99 


»00  <  N^  <    10* 

c        =        144/  o3?6 
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NK  >   10* 

c-d  =  i.  18 

In  addition  to  the  drag  force  the  skin  friction  force  was  calcu- 
lated for  each  layer  and  summed  for  the  particle.   The  friction  force,  a 
result  of  a  minute  boundary  layer  effect  on  the  particle  itself,  was 
relatively  small  but  in  some  cases  significant.   For  spheres  in  low  Rey- 
nolds Number  flows  the  skin  friction  coefficient  can  be  approximated  by: 


/K) 


The  tendency  of  the  velocity  gradient  to  cause  a  moment  on  the  particle 

was  neglected. 

Lift  Due  to  Wall  Proximity. 

Theories  accounting  for  the  effect  of  a  wall  constraint  in  non- 
uniform flow  do  not  exist,  but  this  effect  was  estimated  using  the  exist- 
ing theory  for  uniform  flow  [l2] .   The  effect  of  the  ground  surface  con- 
straint is  to  distort  the  flow  over  the  sphere,  crowding  the  streamlines, 
causing  a  negative  lift  force.   The  theory  is  based  on  kinetic  energy 
considerations,  maintaining  equilibrium  on  a  sphere  in  the  presence  of  a 
wall.   For  a  sphere,  in  a  flow  parallel  to  the  wall,  an  upward  force  is 
required  to  maintain  equilibrium,  indicating  that  the  sphere  must  be 
attracted  to  the  wall. 

The  method  of  calculation  of  this  lift  force  was  similar  to  that  of 
drag.   The  strip  approximation  method  was  used  to  determine  an  average 
local  velocity  over  a  particular  strip.   It  was  then  assumed  that  this 
average  local  velocity  was  that  existing  in  the  vicinity  of  the  stagnation 
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point,  realizing  that  the  stagnation  point  is  shifted  toward  the  region 
of  higher  velocities.   The  lift  then  becomes  a  function  of  the  sphere 
size  and  the  average  local  velocity,  for  a  particle  bedded  on  the  ground 
plane.   The  kinetic  energy  of  a  sphere  in  a  moving  fluid  near  an  in- 
finite fixed  wall  can  be  approximated  by  a  first  order  solution  £.12"]: 

where: 

ro^  =  mass  of  fluid  displaced  by  the  sphere. 
From  Lagrange's  equation  for  equilibrium: 

If  the  sphere  is  to  be  maintained  in  equilibrium,  the  result  is  a  force 
exerted  on  the  sphere: 


fx  =  <//d*  OVbk)  -  s'/i 


A 


F7     =    '/cftrOVdy)     -    ^TA 


7 y  /c7 

Making  the  substitutions,  the  force  exerted  on  the  sphere  in  the  y  or 


vertical  direction  becomes: 


12 


Assuming  that  V,  V  <<*  |. 


where 


we  get: 


^o,-  ^VaTTr* 


^/  =  3//fc^^  "*  ^Vv«- 


/  ~   //fc  \o^  ^k  "    /y 
Assuming  that  the  particle  is  very  close  to  the  ground  such  that; 

y-  r <<  1 

or 

y^r 

the  lift  due  to  the  wall  constraint  becomes: 
For  the  strip  analysis,  summation  yields: 


Lift  Due  to  Velocity  Gradient. 

There  have  been  a  number  of  theoretical  studies  of  the  effects  of 
velocity  gradients  or  shear  on  aerodynamic  forces.  Most  of  these,  how- 
ever, deal  either  with  an  unbounded  flow  or  two  dimensional  flow.   In 
comparing  two  dimensional  solutions  with  three  dimensional  solutions  it 
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was  noted  that  the  difference  is  a  term  in  the  three  dimensional  solu- 
tion describing  the  lateral  component  of  flow.   For  the  axi-symmetrical 
case  dealt  with  here,  this  component  was  assumed  negligible.   A  solution 
for  non-uniform  flow  on  cylinders  by  Murray  and  Mitchell  [^L3~J  was  em- 
ployed in  conjunction  with  a  strip  analysis. 

In  the  strip  analysis,  the  particle  was  divided  into  layers  of  cir- 
cular cylinders  situated  such  that  their  centers  formed  a  line  perpendicu- 
lar to  the  flow.   Beginning  with  the  final  result  of  Murray  and  Mitchell, 
the  dimensionless  shear  velocity  on  each  cylinder  was  found  to  be: 

^'=    -{sme^co^Oi^e)'  y^i    s\^\W   (sme)] 
+  *Aif   C-Dn  §*^>  K,',    CO  cos  <.*.ne> 

4-xLf  C-0n  +  '    X^^C-    IW,   O    sin^+l)9} 
where : 

Kk  CO  =     %<k  l  K2„  C^  ]  ^«  , 


{jo^-  y^  L?  v+-i)  ^K^^^i)]} 
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,f  L^™ 


+  c-01 


and 

fO")  -  Euler's  Constant 
=  0.57721567 

The  dimensional  shear  velocity  is  given  by: 
which  is  of  the  form  (when  squared) : 

Where : 

R*  =   S\nlG  [cos^   C^inG   -    Vrs/  Cosh    ^sia  ©)• 


S  in©  [cosn^S^©')-     '/rv  s^^  Qsin  ©)  ~\ 
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Ui_  is  the  local  velocity  at  the  center  of  each  cylinder,  and  I,  K  and 
K  are  Bessel  functions  of  the  first  and  second  kind,  and  the  derivative 
of  the  second  kind  respectively.   N  is  a  flow  parameter  dependent  upon 
local  velocity,  radius  of  the  cylinder  and  the  vorticity  at  the  stagna- 
tion point,  and  is  defined  as:  N  =  U]~~/(jj     y    •   The  determination  of 
the  flow  parameter  was  approximated  since  the  stagnation  point  on  each 
cylinder  could  not  be  known  in  advance  of  the  calculations.   In  order  to 
keep  the  shear  velocity  in  the  direction  it  was  known  to  flow,  a  lower 
bound  of  N  —  I  could  be  established.   If  <*J0    and  U  are  evaluated  at 
the  geometric  centers  of  the  cylinders,  experimentation  has  indicated  an 
upper  bound  of  N  =  3  0-3]  •   After  analyzing  the  computed  results  of 
shear  velocity  with  values  of  N  ranging  from  1.1  to  10,  a  value  of 
N  =  1.5  was  accepted  as  being  representative  of  the  physical  model. 
Having  found  the  shear  velocity  the  lifting  force  due  to  shear 
effect  was  calculated  for  each  cylinder  and  summed  over  the  particle. 

The  resulting  expression  is: 

( 


-  b^  V*  f    C^-  ^n)s.nG>c/< 
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Where: 


r^  =■   radius  of  cylindrical  strips 
b   "=   width  of  strip 
P\   =   number  of  strips 


Noting  that: 


J        Uj"  Sin©-  c/&        -    O 


We  have : 

Air 


"^  o 

In  evaluating  the  integral  we  have: 

f\lc  s,oe^  -  8M|C-.T^  K',n  CO]  • 
Since  the  last  two  integrals  are  identically  zero.   Further  we  may  write 


Zrr  Jin 


A1  smer/a  *   /     smA  ©•  cosK   CsmeO  r/e 

r7Lrr 
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Since  the  first  and  last  integrals  are  identically  zero. 

Due  to  the  complexity  of  the  above  expressions  and  the  difficulty 
of  integration,  it  was  desired  to  check  this  method  with  another  method 
[.29]  of  a  much  more  approximate  nature  for  calculating  the  lift  due  to 
shear. 

The  alternate  method  of  calculating  lift  assumes  that  the  previously 
discussed  lift  due  to  the  ground  surface  constraint  is  negligible  and 
that  the  lift  forces  on  the  particle  can  be  estimated  by  using  the  solu- 
tion of  Hall  8  for  a  sphere  in  a  two  dimensional  uniform  shear  flow. 
In  addition,  the  approximate  method  assumes  that  the  sphere  rests  on  a 
bed  of  spheres  of  the  same  diameter  and  that  the  pressure  on  the  bed  is 
the  ambient  stream  static  pressure   The  lift  coefficient  then  becomes 
simply  a  function  of  the  local  slope  of  the  velocity  profile,  , 

the  local  velocity  on  the  sphere,  and  the  radius  of  the  sphere.   In  view 
of  the  many  restrictions  placed  on  this  solution  it  was  used  only  as  an 
order  of  magnitude  check  on  the  other  solutions.   Fig.  7  shows  a  compari- 
son of  the  two  solutions. 
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6.    Particle  Entrainment  and  Incipient  Erosion  Rates. 

Only  a  few  of  the  possible  types  of  terrain  to  which  a  mathematical 
analysis  of  the  erosion  rates  might  be  practically  applied,  were  investi- 
gated.  Uniform  particles  ranging  in  size  from  .003  in.  to  .125  in.  in 
diameter  were  investigated.   The  smallest  sizes  represent  loose  dirt 
and  the  larger  diameters  represent  sand  or  sandy  gravel,  with  .125  in. 
being  representative  of  small  gravel  type  rock.   It  was  assumed  that  all 
terrain  was  without  any  moisture  content.   The  densities  of  the  particles 
ranged  from  75  lb/ft3  for  loose  dirt  to  125  lb/ft3  for  gravel.   Table  II 
lists  the  various  terrain  particles  analyzed. 

There  are  three  possibilities  for  entrainment  of  the  ground  particle. 
The  first  is  that  it  may  have  enough  lift  force  on  it  to  be  lifted  directly 
from  the  ground  surface  into  the  free  stream.   The  second  is  that  after 
rolling  motion  is  started  due  to  the  drag  force,  it  may  be  able  to  gain 
the  additional  lift  required  for  entrainment  in  the  free  stream.   The 
third  is  that  with  many  particles  rolling  at  different  velocities  the 
random  collision  of  particles  that  follows  will  bounce  some  of  them  into 
the  inviscid  stream.   Only  incipient  motion  is  within  the  scope  of  this 
work,  hence  it  was  assumed  that  any  particle  which  was  moved  contributed 
to  the  erosion  rate. 

The  net  drag  force  and  the  net  lift  force,  including  both  mechanisms 
of  lift,  were  computed  for  each  of  the  particles,  at  each  of  thirteen 
radial  stations  on  the  ground  surface  for  each  of  the  nozzle  configura- 
tions studied.   To  calculate  the  initial  net  drag  force,  the  previously 
mentioned  ground  surface  model  of  nested  particles  was  assumed  using  a 
coefficient  of  static  friction  of 
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7.   Results  and  Conclusions. 

The  particle  size  criteria  for  lifting  entrainment  and  drag  entrain- 
ment  are  shown  in  Figs.  8,  9,  and  10  in  which  the  particle  size  is  plotted 
versus  a  load  parameter  at  various  stations  with  jet  velocity  and  jet 
height  as  parameters.   The  load  parameter  is  defined  as  the  particle  den- 
sity divided  by  the  dynamic  pressure  in  the  free  stream.   The  lift  en- 
trainment criterion  predicts  that  an  optimum  particle  size,  the  particle 
size  most  prone  to  entrainment,  exists  for  all  radial  stations  and  that 
the  critical  density  increases  with  radial  distance  to  a  certain  ^/[^k,  . 
The  optimum  particle  size  for  lift  entrainment  occurs  because  particles 
smaller  than  the  optimum  are  affected  by  much  smaller  velocities,  and 
particles  larger  than  optimum  are  affected  more  by  a  flow  region  where 
the  shear  gradient,  o  Ly '  ^   approaches  zero. 

The  drag  criterion  indicates  that  extremely  small  particles  of  0.02 
in.  diameter  or  less,  up  to  very  high  densities  will  roll  on  the  ground 
(see  Fig.  8).   As  the  particle  size  is  increased  the  possibility  of 
rolling  decreases  rapidly  up  to  a  particle  diameter  of  0.03  in.  at  which 
size  the  tendency  to  roll  increases  again,  since  the  larger  velocities 
in  the  boundary  layer  begin  to  act  on  the  particle,  until  an  optimum 
diameter  of  approximately  0.06  in.  is  reached.   For  radial  distances 
less  than  WD  =  .1,  velocities  in  the  boundary  layer  are  not  suffi- 
cient to  move  large  particles  of  high  density.   The  range  of  optimum  size 
particle  for  drag  entrainment  of  approximately  .03  to  .09  is  consistent 
for  all  radial  stations.   It  is  of  particular  significance  that  the  range 
of  optimum  particle  diameter  for  both  lift  and  drag  is  very  nearly  equal 
to  the  boundary  layer  thickness  at  the  particular  radial  station  being 
considered. 
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The  incipient  erosion  rates  for  selected  densities,  particle  sizes, 
and  nozzle  configurations  are  shown  in  Figs.  11  and  12.   These  are  found 
by  combining  as  vectors  the  incipient  horizontal  and  vertical  accelera- 
tions of  the  particles.   This  was  assumed  to  be  proportional  to  the 
incipient  erosion  rate.   The  erosion  rate,  normalized  by  the  maximum 
rate,  is  plotted  versus  the  nondimensional  radial  distance.   For  all  con- 
figurations and  particles  the  erosion  profile  is  maximum  in  the  area 
.5  *C  ^/t-s  <^.0.  The  maximum  erosion  for  larger  particles  occurs  near 
^/xa   —  1,  and  then  drops  off  sharply  with  increased  radial  distance. 
As  the  particles  become  smaller  a  maximum  erosion  occurs  closer  to  the 
perimeter  of  the  jet  nozzle  and  a  second  maximum  begins  to  occur  at 

/^N  ^^  1.3.   This  occurs  after  the  transition  region  has  begun  to 
form  and  where  the  onset  of  turbulence  can  affect  the  particles.   Fig.  13 
shows  the  incipient  rate  of  lift  entrainment  by  itself  and  indicates  that, 
in  general,  the  shape  of  the  erosion  profile  can  be  attributed  to  the 
lift  forces. 

The  results  shown  in  the  graphs  represent  only  a  few  configurations 
selected  from  the  large  amount  of  data  collected.   In  viewing  all  of  the 
results  and  comparing  the  many  combinations  of  configurations,  it  was 
found  that  the  variation  of  ground  erosion  rates  was  principally  a  func- 
tion of  the  dimensionless  radial  distance  from  the  jet  centerline  and 
the  size  of  the  particle  for  a  given  nozzle  height,   /n^r^  •   The  ero- 
sion rate  was  nearly  proportional  to  the  particle  density  except  near  the 
stagnation  point,  where  the  erosion  rates  were  very  small. 

From  these  results  it  may  be  concluded  that  for  ground  particles  of 
less  than  .04  in.  diameter  the  incipient  erosion  rate  has  two  maximum 
points.   One  in  the  vicinity  of  the  perimeter  of  the  jet  nozzle,  and  the 
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second  at  a  radial  distance  approximately  two  diameters  from  the  stagna- 
tion point.   For  very  small  particles  and  low  jet  velocities  the  second 
maximum  is  the  most  significant.   In  all  other  cases  where  two  maximum 
points  exist,  the  first  is  the  most  significant.   For  ground  particles 
whose  diameters  are  greater  than  .04  in.  a  single  maximum  point  exists 
near  the  perimeter  of  the  jet  nozzle.   As  the  particle  size  increases 
the  influence  of  the  turbulent  region  decreases. 

It  was  concluded  that  particles  of  approximately  .02  in.  diameter 
or  less  would  entrain  due  to  drag  even  at  very  high  densities  or  very 
low  velocities.   For  particles  larger  than  approximately  .02  in.,  an 
optimum  size  of  approximately  .06  in.  exists  up  to  radial  distances  of 
four  nozzle  diameters  for  drag  entrainment.   An  optimum  particle  size 
most  susceptible  to  entrainment  due  to  lift  exists  in  the  range  of  .07 
to  .12  in.  diameter  for  radial  distances  up  to  four  nozzle  diameters. 

Also,  it  may  be  concluded  that  varying  the  nozzle  diameter  had  very 
little  effect  on  the  dimensionless  erosion  rate;  hence  in  the  final  analy- 
sis the  nozzle  diameter  was  not  considered  as  a  parameter.   Moreover, 
while  the  actual  erosion  rates  were  approximately  linearly  dependent  upon 
the  jet  velocity,  the  dimensionless  erosion  profile  appeared  to  be  inde- 
pendent of  jet  velocity,  with  the  exception  of  the  order  of  importance 
of  the  maximum  points. 
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VELOCITY  PROFILES  FOR  LAMINAR  AND  TRANSITION 
BOUNDARY  LAYER,   Z/DN  =r  .5 
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APPENDIX 

FORTRAN  PROGRAMS 

LIFT  FORCES: 

Input. 

Program  is  generated  from  the  polynomials  describing 

the  velocity  profiles  of  the  boundary  layer,  for 

various  radial  stations, 

MM  Number  of  terms  of  the  polynomial 

PCOZ  Coefficients  of  the  polynomial 

QSH  Nondimensional  shear  velocity 

VJ  Jet  velocity 

DN  Nozzle  diameter 

RG  Radial  distance  from  jet  centerline 

ZG  Nozzle  height  above  ground 

YHVX  Distance  to  have  maximum  velocity  in  turbulent 

boundary  layer 

VMX  Maximum  velocity  in  turbulent  boundary  layer 


Intermediate  Results. and  Output 

DP  Particle  diameter 

VI  Inviscid  velocity 

YND  Nondimensional  distance  above  ground 

XND  Nondimensional  roots  of  velocity  profile 

polynomial 

DNS  Particle  density 

DVX  Derivative  of  polynomial 

SDVX  Slope  of  velocity  profile  •- 

DUDY  Slope  of  velocity  profiie  - 

FLSC  Shear  lift  {cylinder  method) 

SLFT  Shear  lift  (sphere  method) 

PLFT  Lift  due  to  wall  constraint 

PLAC  Acceleration  of  particle  in  vertical  direction 

Subroutine  Root  finds  the  roots  of  the  velocity  profile 
polynomial 
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611 
612 
400 
401 

4o3 

303 

302 
661 

640 
61 
19 
60 

18 

20 

22 

64i 
644 
645 

647 

643 
654 
655 


SUBROUTINE  LIFT(     MM,PCOZ,    QSH , V J , ON , RG, ZG » NCC , YHVX , VMX  ) 

DIMENSION  DP(10),RP<10),        YGU ( 4 0 ) , YGS [ 40 ) , YL ( 70 ) »  YG ( 7 0  )  » 

AR(7o),VL(7o)»PCOZ(3o),XL(4o>,FX(3o)»NlY(5),DVX(30>» 
!  DGb(5),PCOX(10),R(100) > YND ( 12  0 ) , XND ( 120 ) , YNDDC50)  , XNDD<50># 
I  VLD(5o ) ,  YGDD(2o ) 

UP(D  =  .06 

DP(2)= . 04 

DP(3)=.Q2 

UP(4)=. 01 

DP<5)= .007 

DP(6)=.125 

DP(7)  =  .003 

RGN=RG/DN 

ZGN=ZG/DN 

PRIN1  611,  DN 

FORMAT*/// 

PRINT  612, 


VJ 


22HN0ZZLE  DIAMETER  (FT)  =   F5.2  ) 

FORMAT(  26HN0ZZLE  VELOCITY  (FT/SEC)  =   F7.1  ) 

PRINT  400*RGN 

F0RMAT(26HN0NDIMEN  RADIAL  DISTANCE  =   F20.10) 

PRINT  401,  ZGN 

F0RMAT(26HN0NDIMENJ  HEIGHTH  OF  NOZZLE  =   F20.10) 

DO  115  I  1=1,7 

PRINT  403,  DP  < I  I  ) 

F0RMAT(//23HPARTICLE  D  I  AMETER ( I N . ) =  F2Q.10>- 

RP( I  I )=DF( I  I  )/2. 

GNU=3.745E-7/2.3769E-3 

IF(NCC-25)3o2#303,302 

VI=VMX 

GO  TO  70  0 

I  F(ZG/DN-. 5)64Q,640,661 

KZG=ZG/DN 

GO  TO(640,64l,640,643) ,KZG 

IF(RGN-.35)60,60,61 

IFCRGN-. 75)18,18,19 

IF(RGN-1.2)20,20,22 

VI=VJ*S0RTF(ABSF(1.88*RGN**2-.225*RGN)) 

GO  TO  700 
VI=VJ*SQkTF(-3,2*RGN**2+5,13*rgn-1.26) 

GO  TO  700 
VI=VJ*SQRTF(-2,68*RGN**2+5.48*RGN-1.82) 

GO  TO  700 
VI=VJ*SQRTF(2.25*RGN**2-9.45*RGN+8.96) 

GO  TO  700 

IF(RGN- .35)6q,6o ,644 

IF(RGN-1. 0)645,645,647 

VI=VJ*SQKTF(-.Q985*RGN**2+1.13*RGN-.24D 

GO  TO  700 
VI=VJ*SQRTF(.077*RGN**2-.641*RGN*1.354) 

GO  TO  700 

IFCRGN-. 40)62, 62, 654 

IF(RGN-1. 2)655, 655, 657 

Vl=VJ*SGRTF(-.5*RGN**2+1.2*RGN-.3) 


92 


657 

62 
700 
016 


1 
100 
901 
705 

704 

?00 

L04 

504 
500 

511 
501 


502 
506 


507 

508 

71 

105 

573 

32 

31 
706 

44 

45 
43 


40 


GO  TO  7 
V  I  =VJ*S 
GO  TO  7 
VI=VJ*S 

ANCL=10 
TH=DP(  I 
NN=ANCL 
N  A  =  (  N  N  + 
R(D=RP 
DO  1  K  = 
AA=K-1 
XLC=AA* 
R(K)=SU 
IF(NCC- 
IFCNCC- 
YND(  1) 
GO  TO  1 
YND(  1) 
GO  TO  1 
AAA=VI/ 
YNDC  1) 
XND(  1  )  = 
IF(NCC- 
IF( YND( 
MM  =  8 
DO  511 
PCQX(  11 
GO  TO  5 
MM  =  4 
DO  5q2 
N=I+8 
PCOX(  I  ) 
DO  507 
EX1=MM- 
F  X  (  K  W  )  = 
XND(  1) 
GO  TO  6 
DO  71  K 
EX1=MM 
FX(KZ)= 
XND(D  = 
GO  TO  6 
CALL  HO 
IF(XND( 
XND(D  = 
VLll)=V 
A  R  T  =  0  • 
DO  43  J 
IF( JZ-1 
AKt  =  2-  * 
GO  TO  4 
ARE  =  4.* 
A  R  T  =  A  R  T 
SPA=RP( 
ADS=SPA 
ADS=ADS 
ACOR=l. 
FLS=0- 
DO  39  I 
IFCIZ-1 
FL=-TH/ 
GO  TO  3 


00 

Q R T F  (  .(j404*RGN**2-.3l8*RGN-»-.744} 

00 

Q R T F  (  .625*RGN**2-»-.o23*RGN) 

1  . 

I  )/ANCL 

i)/2 

(II) 

2,NA 

TH 

R  IF(KP( I  I  )**2-XLC**2> 

26)901.900,900 

25)704*705,900 

=  R  (  1J_/ JTH V X 

04 

=(R(  1)/<12.*DN))*SQMTF(VJ*DN/GNU) 
04 
RG 

=R(1)*SQRTF(2.*AAA/GNU)/12. 
0.0 

25)508.504,508 
1)-. 15)500  ,500,501 

I  I  =  1 ,  M  M 
)=PCOZ( IZ) 
08 

1=1, MM 

• 

=PCOZ(N) 
K  V  = 1 ,  M  M 

PCOX(KW)*VND(  1)**EX1 

=XND(  1)+FX(KW) 

73 

Z=1,MM 

-KZ  +  1 

PCOZ(KZ)*YNJD(D**tXl 

XND(1)+FX(KZ) 

73 

Oi (1,MM,PCOZ,XND.YND.NCO 

D-l.  0  )31,31,32 

.9999 

I*XM)(1) 

Z=1,NA 

)44. 44, 45 

R( JZ)*3.14159*TH 

3 

R( JZ)*3.14159*TH 

+  ARb 

I  I  )**2*4,*3. 14159 

-AR1 

/SPA 

+  ADS 

Z=1,NA 

)  40.40,41 

2.*2.3769E-3*R(IZ)*VL(1)**2*QSH/144. 

9 
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LIFT  OF  PARTICLE  BY  SPHERE  METHOD 


*AAA/GNU)/12 


39  FLS=FLS+FL 

FLSC=FLS*ACOR 
APPROXIMATE  SHEAR 
130  YGDDd  )=R(D 

IF(NCC-25)132,77,129 
11     YNOD(l)=YGDD(l)/YHVX 
GO  TO  133 
129  YNDD( 1)=YGDD(1)*SQRTF(2 
GO  TO  133 

132  YNDD(1)  =  <YGL)D(1)/<12.*DN))*SQRTF(VJ*DN/GNU) 

133  SDVX=0.0 

DO  74  KD=1»MM 
MX2=MM-KD 
EX1=MM-KD+1 
IF<NCC-25)75, 79,75 

79  PCOZ(KD)=PCOX (KD) 

75  IF(KD-MM)lB0fl61#lti0 

180  DVX(KD)=EX1*PC0Z(KD)*YNDD(1)**MX2 
GO  TO  14 

181  DVX(KD)=PCOZ(KD) 
74  SDVX=SUVX*DVX(KD) 

IF(S0VX)182»l83,183 

182  1F(NCC-25)184#183,184 
184  SDVX=0.0 

IF(NCC-25)190,183,192 

183  DUDY=SDV**VMX/YHVX*12. 
GO  TO  193 

192  DUDY=SDVX*SQRTF(2.*AAA/GNU)*VI 
GO  TO  193 

190  DUDY=SDVX*SORTF(VJ*DN/GNU)*VI/DN 

193  ZON=     K(1)/VL(D*DUDY/12, 
CLF=0.345*0.557*ZON+0.879*ZQN**2 
SLFT=CLF*2.37  69E-3*VL(1)**2/2.*3.14159*R(1)**2/144. 

PRINT  800 

800  FORMAT(14HI\ONDIMEN  SHEAR 
1  SHEAR, 5X,17HLIFT  DUE  TO 
2FILE  SLOPE) 

PRINT  801 

801  F0RMAT(3X,bHVEL0CITY,HX 
1    8H(SPHEKE)/) 

PRINT  B02,GSH,XND(1),FLSC,SLFT,ADS,SDVX 

802  FORMAT(6E20  .10) 
GO  TO  52 

5  PRINT  200 
200  F0RMAT(/5HEKR0R/) 
52  NNI=59 

UN=NNT 
603  QGB(1 )=0 . 

YH=0. 

YLM=NNT-1 

YINC  =  DP(  1  I  )/YLM 

DO  169  J=1,NNT 

IF( J-D666, 111,112 

111  YG( J)=0 . 
GO  TO  169 

112  YG( J)=YH+YINC 
YH=YG( J) 

169  CONTINUE 

DO  113  K=2,NNT 
YND(D  =  0. 
IF(NCC-26)116,902,902 


5X,17HN0NDIMEN  VELOC I T Y , 5X , 17HL I  FT  DUE  TO 
SHEAR, 5X.19HAREA  DIFFERENCE  0/0 , 3x , 13hPR0 


11HON  C YL I NDER»12X,10H( CYLINDER) »13X, 


9k 


902 

116 
301 

300 
113 
671 


720 
721 

811 

701 

702 

707 

708 


37 
672 


33 
38 

51 
35 
34 
95 

FI 
96 


114 


731 
732 


YND(K) 
GO  TO 
IF(NCC 
YND(K) 
GO  TO 
YNIMK) 
CONT  IN 
DO  672 
X  N  D  (  N  X 
IF(NCC 
IF( YND 
MM  =  8 
DO  811 
PCOX(  I 
GO  TO 
MM  =  4 

DO  702 
N=  1  +8 

PCOX(  I 

DO  707 

EX1=MM 

FX(KW) 

XND(NX 

GO  TO 

DO  37 

EX1=MM 

FX(KL) 

XND(NX 

CONT IN 

DO  95 

IF(XND 

IFCKC- 

XND(KC 

GO  TO 

X  N  D  (  K  C 

IF(XND 

X  ,M  D  (  K  C 

CONTIN 

ND  AVE 

SUMV=0 

DO  114 

VL(L)= 

SUMV=S 

VLAV=S 

PLIFT= 

PL  IF  I  = 

V0L=4. 

DO  730 

GO  TO( 

DNS=75 

GO  TO 

DNS=90 

GO  TO 


=  YG(»\)*SURrF(2.*AAA/GNll>/12. 

113 

-25)300,301,902 

=YG(K)/YHVX 

113 

=(YG(K)/(12.*DN))*SQRTF(VJ*DN/GNU) 
UP. 

N  X  =  1  ,  N  N  T 
)  =  0. 

-25)708,720,708 
(NX)-. 15)721  ,721,701 

I  Z  =  1 ,  M  M 
Z)=PCOZ(  IZ) 
708 

1=1, MM 

)=PCOZ(N) 

Kw=l,MM 
-KW 

=PCOX(KW)*YND(NX)**EXl 
)=XND(NX)*FX(KW) 
672 

K  L  =  1 ,  M  M 
-KL  +  1 

=PCOZ(KL)*YND(NX)**EXl 
)=XND(NX)+FX(KL) 
UE 

KC=1#NNT 
(KC) )33,35,35 
1)38,38,51 
)  =  0.0 
35 

>=XNU<KC-1  ) 
(KC)-i. 0)95,34,34 
)=.9999 
UE 
RAGE  VELOCITY 

• 

L=1,NNT 
V  I  *  X  N IJ  (  L  ) 
UMV-fVL(L) 
U  M  V  /  U  M 

(2.3/69E-3/2.)*VLAV**2*3.14l67*(3./8,)*RP(II)**2/144l 
-PL  I  I-  1 
/3.*3.14159*KP(  I  I  )**3 

M  U  =  1 ,  4 
731, 732, 733, 734), MO 

735 

735 


95 


733 

734 
735 

736 


842 
843 

110 


120 


840 
730 

666 
667 

115 
668 


199 
201 
106 


101 

102 
103 

100 


104 


DNS=105. 
GO  TO  73 
DNS=125. 
PW I =DNS/ 
PRINT  73 
FORMAT  (/ 
PL  I  OT  = 
PMAS  =  P*IT 
PLAC=PLT 
IF(PLAC) 
PLAC=0 . 0 
PRINT  11 
FORMATC/ 
.  26HNET 
PRINT  12 
FORMAT(7 
.AL  DIREC 
PRINT  64 
FORMAT(F 
CONTINUE 
GO  TO  11 
PRINT  66 
FORMAT(/ 
GO  TO  66 
CON1 INUE 
REIURN 
END 

SUBROUTI 
DIMENSIO 
DO  210  L 
IF(NCC-2 

IF(YND(L 

XND(D=0 

SFX=0. 

SDVX  =  0  . 

DOl00N=l 

MX1=MM-M 

MX2=MM-N 

EX1=MX1 

FX(N)=PC 

IF(N-MM) 

DVX(N)=E 

GO  TO  10 

DVX (N)=P 

SFX=SFX+ 

SDVX=SDV 

RMX=(SFX 

IF(ABSF( 

XND(D=X 

MLM=1000 

IF(XND(L 

GO  TO  21 


1728. *VOL 

6,  DNS 

1BHPARTICLE  DENSITY  =  F6.1) 

PLIFT-PWT+FLSC 
/32.16 
OT/PMAS 
642*843,843 


0 


' 


22HAVERAGE  LOCAL  VELOC I T Y , 3X , 16HL I F T  ON  PARTICLE, 3X, 

LIFT  FORCE  ON  PART  I  CLE , «X , 2lHP ART  I  CLE  ACCELERATION) 

0 

X,6H(FT/SEO ,7X,22HUUE    TO    WALL    CONSTRA I  NT , 30 X , 2lH I N    VERTIC 

I  ION/) 

0/ VLAV,PLIFT,PLTOT,PLAC 

20.10>E2Q.10,5X,F2Q.10,5X,F2Q.10/) 

5 
7 

5HERR0R) 
8  • 


NE    ROOT(NNT,MM,PCOZ, YND.XND,NCO 

N  YND(70)»XND(70),PCOZ(30)»DVX(30)#FX(30) 

=1,NNT 

4)199,199,201 

)-. 425)200, 201, 201 

•  1 


,  MM 
+  1 


OZ( 
101 
XI* 
3 

COZ 
FX( 
X  +  D 
-YN 
RMX 
in*D  ( 
/MM 

)-2 

0 


N)*XND(L)**MX1 

,102,101 

PC0Z(N)*XND(L)**MX2 

(N) 

N) 

VX(N) 

D(L) )/SDVX 

)-,01     )108,108,104 

D-RMX 

,**MLM)106,204,204 
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108  IF(NCC-^>109, 109,210 

109  IF(NNT-?)210,210,11Q 

110  IF(XND(L)-1. 0)^10,203,203 

203  IF(YNU(L)-2. 0)204, 204. 205 
205  XND(D  =  .9999 

GO  TO  210 

204  IF<L-l>3o4,304,3o5 

304  XNO(D=0  .  0 
GO  TO  210 

305  XND(D=XNDIL-1) 
GO  TO  210 

200  GO  TO  210 
210  CONTINUE 

RETURN 

END 

END 


DRAG  FORCE: 

Input. 

Program  is  generated  in  the  same  manner  as  the  lift 
program. 

Intermediate  Results  and  Output. 


DELT 

XDEL 

FDT 

FFDT 

PWT 

PDAC 


Boundary  layer  thickness 

at  top  of  boundary  layer 
Drag  due  to  pressure  force 
Drag  due  to  friciton 
Particle  weight 
Horizontal  acceleration  of  particle 


SUBROUTINE  DRAG(PCOZ,MM,NCC,VJ,RG,ZG.D<N*YHVX,VMX) 

DIMENSION  DP(10),UN(40),RP(10),YGU(40),YGB(40),YL(70),YG(70>* 

1  AR(70)*YND(70).XND(70),VL(70),PCOZ(30),XL(40),FX(30), 

2  XDELC5) , YMAX(5),PCOX(10  )  , ARF(70  ) 
DPU)  =  .06 

DP(2)=,04 
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DP(3)=.02 
DP(4)=,01 
DP(5)=.007 
DP(6)=.l25 

DP(7)=,003 
N  N  =  3 1 

RGN=RG/DN 
ZGN=ZG/DN 
PRINT  611. 

611  F0RMAK22HN 
PRINT  612, 

612  F0RMAT(26HN 
PRINT  4oo»R 

400  F0RHAT(26HN 
PRINT  401, 

4oi  F0RMAT(28HN 
DO  115  I  1=1 
PRINT  403. 

4o3  FORMAT*  23h 
RP( I  I )=DP(  I 
UNI=NN 
TK=RP<  I  I  )/( 
YLC=0. 
SUMAR=0 ♦ 
NAl =NN-1 
DO  3  N=1,NA 
YL(N)=YLC+T 
AM  =  N 

XL(N)=SQRTF 
AR(N)=2.*XL 
ARF(N)=^.*X 
IF(N-1)8#5# 

5  SUMAR=SUMAR 
GO  TO  6 

8  SUMAR=SUMAR 

6  YGU(N)=RP(  I 
YGb(N)=RP(  1 

3  YLC=AM*TK 
GNU=. 000000 
IF(NCC-25)3 
303  VI=VMX 

GO  TO  700 
302  IF(ZG/DN-.5 
661  KZG=ZG/DN 

GO  TO(640,6 

640  IFCRGN-.35) 
61  IFCRGN-.75) 

19  IF(RGN-1.2) 
60  VI=VJ*SQRTF 

GO  TO  700 
18  VI=VJ*SQRTF 
GO  TO  700 

20  VI=VJ*SQRTF 
GO  TO  700 

22  VI=VJ*SORTF 

GO  TO  700 

641  lFlRGN-,35) 

644  '  IFCRGN-1. 0  ) 

645  VI=VJ*SQRTF 
GO  TO  70  0 

647  VI=VJ*SQRTF 
GO  TO  700 


UN 

OZZLE  DIAMETER  (FT)  =   F5.2) 

VJ 

OZZLE  VELOCITY  (FT/SEC)  =   F7.1  ) 

GN 

ONUlMEN  RADIAL  DISTANCE  =   F20.10) 

ZGN 

ONDIMEN  HEIGHTH  OF  NOZZLE  =   F2Q.10) 

.7 

DP  (  I  I  ) 

PARTICLE  DIAMETER-IN.  =   F2Q.10  ) 
I  )/2. 

UNl-,5) 


T   • 
K/2. 

(RP( I  I )**2-YL(N)**2) 

(N)*TK 

L(N)*TK 

8 

+3.1416/2. *ARF(N> 

+3.1416*ARF(N) 
I  )+YLC 
I  )-YLC 

3745/. 0023769 
02.303.302 


)640.640.661 

41.643.643)  .KZG 
60 ,60 ,61 

18,18.19 
20.20.22 
(RGN**2+.0175*RGN) 

(-3,2*RGN**2+5.13*RGN-1.26) 

(-2,68*RGN**2t-5.48*RGN-1.82) 

(.256*RGN**2-1.48*RGN*2.3ll) 

60.60, 644 

645.645,647 
(-.0985*RGN**2+1.13*RGN-.24D 

(.077*RGN**2-.641*RGN*1.354) 


98 


643 
654 
655 

657 

62 
700 


24 
26 

27 

28 
977 

900 


99 
301 

978 

29 

672 


804 
800 

811 

801 

802 
806 

807 

808 


809 

805 

75 

33 
72 

73 

670 

34 

32 

74 
79 

76 


IFCRG 
IF(RG 
VI=VJ 
GO  TO 
VI  =  VJ 
GO  TO 
VI=VJ 

N  N  T  =  2 
D029N 
I  F  (  N  - 
YGCN) 
GO  TO 
YGCN) 
I  F  (  N- 
ARCN) 
GO  TO 
AR(N) 
IFCNC 
AAA  =  V 
YND(N 
GO  TO 
IF(NC 
YND(N 
GO  TO 
YNDCN 
CONTI 
DO  80 
XND(N 
IFCNC 
IF(  YN 
MM  =  8 
DO  81 
PCOX( 
GO  TO 
MM  =  4 
DO  80 
PCOX( 
DO  80 
EX1  =  M 
F  X  (  K  W 
XNDCN 
GO  TO 
DO  80 
EX1  =  M 
FX(KY 
XNDCN 
CONTI 
DO  32 
IFCXN 
1FCKC 
X  N  D  C  K 
GO  TO 
X  N  D  (  K 
IFCXN 
XNDCK 
CONTI 
IFC  I  I 
IFCNC 
DELT  = 
GO  TO 
XDEL  = 
CALL 


N-.40) 
N-1.2) 
♦SQRTF 

700 
*SQRTF 

70  0 
*SQRTF 
*NN-3 
=1,NNT 
NN)24, 
=  Y  G  B  C  N 

27 
=  Y  G  U  (  N 
NN)27, 
=ARCNN 

97  7 
=ARCN- 
C-26)9 
I/RG 
)=YG(N 

29 
C-25)9 
)=YG(N 

29 
)=( YGC 
NUE 

5  NW  =  1 
W)  =  0.0 
C-25)6 
D(NW)- 

1  IZ  =  i 
IZ)=PC 

808 

2  1=1, 

I )=PCO 
7  KW  =  1 
M-KW 
)=PCOX 
W) =XND 

805 
9  KYsi 

M  -  K  Y  ♦  1 
)=PCOZ 
W )=XNU 
NUE 

KC  =  1, 
DCKC)  ) 
-1)72, 
C)=0.0 

670 
C)=XND 
DCKC)- 
C)=.99 
NUE 

-1)74, 
C-25)7 
4.6*12 

80 
0.99 
ROOTC1 


62,62, 654 

655, 65d, 657 

C-.5*RGN**2+1.2*RGN-.3) 

(,04o4*RGN**2-.3l8*RGN+.744) 
C.625*RCiN**2*.o25*RGN) 

26,26 
N-N) 

-NN+2) 

27,28 

-N) 

NN  +  2  ) 
9,900,900 

)*SQRrFC2.*AAA/GNU)/l2. 

78,301,900 
)/YHVX 

N)/Cl2,#DN))*SQRTFCVJ*DN/GNU) 

,NNT 

08.804,808 
.15)800,800,801 

,MM 
OZ(IZ) 


MM 

ZC 1+8) 

,MM 


CKW)*YND(NW)**EX1 
(NW)+FXCKW) 

,MM 

<KY)*YNU(NW)**EXl 
CNW)+FXCKY) 

NNT 

33,670,670 

72,73 


(KC-1) 
1.0)32,32,34 

99 

74,80 

8,76,79 

./SQRTFC2.+AAA/GNU) 


,MM,PCOX,XDEL,YMAX,NCC) 


99 


DELT  =  YMAX(D*YHVX 

GO    TO    80 
78    XDbL(l)=.99 

CALL    R  0  0  I  ( 1 , MM  #  PCOZ  *  X  DEL  *  Y  M  AX  #  NCC ) 

DELT=YMAX(1)*DN/SQRTF(VJ*DN/GNU)*12. 

80  GO  TO  31 
31  D035K=l,NNi 

35  VL(K)=VI*XiMD(K) 
82  NMID=(NN !  «■  1  )  /  2 

RN  =  VL(NM1IJ)*DP(I1)/(12.*GNU) 
IFCRN-100. )730,730, 731 

731  1F(RN-10000.  )732,732,733 
730  CD=30 ./(RN** .699) 

GO  TO  734 

732  CD=1 .44/(RN**.o396) 

GO  TO  734 

733  CD=1.18 

734  FDT=0. 
D036K=1,NNT 

FD=2.376  9E-3*CD*VL(K)«*2*AR(K)/288. 

36  FDI=FDT+FD 

ARS  =  4.*3.14159*RP( I  I  )**2 

ARD=ARS-SUMAR 

ARD=ARU/AHS 

ACOR=l.+ARD 

CF=     .664/SURTF(RM) 

FFDT=0. 

DO40M=l,NNl 

FFD=2.3769b-3*CF*VL(M)**2*ARF(M)/288. 

40  FFDT=FFD(+FFD 
CFB=4.*CF 

TOL=ABSF(CKB*. 00001) 
CALL  CUBE(CFB,T0L,CFB3) 
CBU=0.1/CFB3 

FBDT=0 . 

UO  42  MD=1*NNT 

FBD=2.3769b-3*CBD*VL(MU)**2*AR(MD)/288. 
42  FBDT=FBDT+FBD 

PRINT  38 
38  F0RMAT(/4X, 16HPRESSURE  DRAG-LB, 4X » 16HFR I CT I  ON  DRAG-LB , 8X , 12HBASE  I 
1RAG-LB,2X,24HB0UNDARY  LAYER  TH I CKNESS , 2X , 31HSPHERE  APPROXIMATION  I 
2IFFERENCE//) 
PRINT  41,  KUT,FFDT,FBDT,DELT,ARD 

41  FORMAT(   3E20.10/2F20.10) 
V0L=4./3.*3.14159*RP( II )**3 
DO  720  MO=l*4 

GO  TO(72l,722,723,736) ,MO 

721  DNS=75. 
GO  TO  724 

722  DNS=90. 
GO  TO  724 

723  DNS=105. 

GO  TO  724 
736  DNS=125. 

724  Pwl=DNS/1728.*VOL 
PRINT  725, DNS 

725  F0RMAT(/18HPARTICLE  DENSITY  =   F6.D 
TOIL=(FDT+FFDT-,7  07*PWT)*ACOR 
PRINT  620,  TOTL 

620  F0RMAT(/34HNET  FORCE  IN  DIRECTION 
PMAS  =  PWT732.l6 
PDAC=TOTL/PMAS 

100 


OF  MOTION  =   E20.10) 


- 


IF(PDAC)842/843,843 
PDAC=0 . 0 
PRINT  621,  PDAC 
F0KMAT(/56HPARTICL£ 
1=   F20.10) 
720  CONTINUE 
CONTINUE 

IF(NCC-24>668»668,5>99 
CONTINUE 
RE1URN 
END 


842 
843 

621 


115 

599 
668 


ACCELERATION  IN  HORIZONTAL  DIRECTION  (FT/SEC) 


101 


